Copyright © 1979 American Telephone and Telegraph Company 

The Bell System Technical Journal 

Vol. 58, No. 2, February 1979 

Printed in U.SA. 



Properties of Kruithof s Projection Method 
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(Manuscript received July 24, 1978) 

In 1937, J. Kruithof introduced a scheme for projecting from mea- 
sured point-to-point teletraffic data to some future values, based upon 
estimates of total originating and terminating traffic only. This study 
seeks to give a unified picture of Kruithof 's projection method and its 
generalizations, with some practical details and recommendations 
for implementation. The main text deals with existence and conver- 
gence testing, treatment of ill-conditioned or slowly convergent cases, 
and various extensions of the basic method. An appendix includes 
proofs of existence, uniqueness, convergence, and continuity. 

I. INTRODUCTION 

J. Kruithof s method 1 for projecting from measured point-to-point 
teletraffic data qy to some future values py is based upon estimates of 
total originating and terminating traffic only. While the original pub- 
lication (in Flemish) did not receive as much attention as it may 
deserve, the idea was good enough that it has been independently 
reinvented numerous times in the intervening years. Related tech- 
niques have turned up in economics, statistics, biophysics, pattern 
recognition, and vehicular traffic studies, for instance. Such repetition 
largely seems due to a scientific "Babel" effect: workers in different 
technical disciplines can no longer read each other's work and recognize 
the same problem in a new context. 

While Kruithof showed that his method had certain properties 
which are clearly desirable in a projection scheme, he did not investi- 
gate the underlying mathematical problems. Subsequent workers, such 
as Bear, 2 Kullback, 3 Sinkhorn, 4 ' 5 Theil, 6 and particularly Csiszar, 7 have 
thrown much light on these matters. This paper seeks to give a unified 
picture of Kruithof 's method and its many generalizations, with some 
practical details and recommendations for implementation. Much of 
the more intricate mathematics is relegated to an appendix, including 
proofs of existence, uniqueness, convergence, and continuity. These 
are cited as needed in the main text, which contains information on 
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existence and convergence testing, treatment of ill-conditioned and 
slowly convergent cases, and various extensions of the basic method. 
The final comments propound a rationale for Kruithof projection in 
the context of Bell System planning. 

II. KRUITHOFS BASIC METHOD 

In 1937, J. Kruithof 1 proposed a technique for predicting the point- 
to-point traffic pij in a given year from a number of originating points 
i = 1, 2 • • • M to a number of terminating points / ■ 1, 2 • • • N. The 
units could be calls, trunks, or erlangs, for instance, as long as it makes 
sense to add up various entries in the matrix p = [py]. It is assumed 
that the corresponding traffic matrix q = [qij ] is known for some other 
year, while total traffic 6, at each point i and dj at each j have already 
been estimated by some external means to yield: 

bi = Y i Pij (1) 

j 

dj=lPij. (2) 

i 

Then Kruithof 's formula for projecting p from q is: 

p i j=q ij E i Fj. (3) 

The "growth factors" Ei and Fj for the originating and terminating 
points are implicitly defined, and must be computed by solving (l)-(3) 
simultaneously. 

Kruithof recommended that (l)-(3) be solved by starting with the 
estimate p = q, then alternately normalizing the rows of p to satisfy 
(1) and the columns to satisfy (2), until it stops changing. In practice, 
this scheme suffers from a tendency to accumulate roundoff error. A 
mathematically equivalent procedure with better numerical properties 
would be to substitute (3) into (1) and (2), solving for Ei and Fj to 
obtain: 

Ei = &,•/£ qijFj (4) 

j 

F j = dj/lq ij E i . (5) 

i 

Starting from an arbitrary estimate, such as Fj ■■- 1 for ally, (4) and (5) 
may be evaluated alternately until p converges. 

Various questions arise naturally in connection with this projection 
method: 

(i) Under what conditions do (l)-(3) possess a solution p *? 

(«*) Can there be more than one solution for p * ? 

(Hi) How does p * vary with the estimates of total traffic? 

(iv) Does the iteration converge, and if so, to what? 

518 THE BELL SYSTEM TECHNICAL JOURNAL, FEBRUARY 1 979 



(v) When should we stop iterating a given case? 
(vi) Are there valid generalizations of this scheme? 

Many of these points are treated at considerable length in the appen- 
dix, which discusses a generalized problem: find a probability distri- 
bution P which satisfies arbitrary linear constraints, such as (l)-(2), 
and is related to a given distribution Q by a product formula, such as 
(3). To make the connection in the present case, our first step is to 
divide p by the sum p of all its elements, reducing it to a joint 
probability P ■ p/p that a call, trunk, or other increment of traffic is 
from i toy. Now (1) and (2) yield relations: 

B, = 6,/p = £ P„ d') 

j 

D J ^dj/p = lP iJ , (2') 

i 

with Bi and Dj the marginal distributions of traffic on i and/ Similarly, 
q is divided by the sum q of its elements to get the joint distribution 
Q = q/q. Let the events {e) in the appendix be the set of pairs e = 
(i,j) and the constraints {c} corresponding to (29) be (l')-(2') for all 
points i and/ Then taking Et = Vip/q and Fj = Wj puts the projection 
formula (3) in exactly the product form (32): 

P u = QijEiFjq/p = QijViWj, (3') 

so that we have a case of the general Kruithof problem defined in the 
appendix. 

Now the existence conditions from the appendix show that there is 
a solution p* of the form (3) if and only if (1) and (2) have some 
solution pij that vanishes for each q$ that vanishes and is positive 
whenever qy is positive. The uniqueness results show there is at most 
one solution p*. The solution is continuous in all 6, and dj whenever it 
exists. The iteration on (4) and (5) can be recognized as an example of 
a relaxation procedure, as discussed and analyzed in the appendix. If 
(l)-(2) possess a solution py that is zero for each qy that is zero, the 
iteration will converge to some p with these same properties. The limit 
may not be of the form (3) though, since p,y can also vanish for some 
qij > 0. When a solution p* to (l)-(3) exists, however, the iteration 
converges to it uniquely. This explains the resistance of (4)-(5) to 
roundoff effects, since such perturbations die out in the process of 
converging. In the special case that q tJ is symmetric and 6, = d, for all 
i, uniqueness shows thatp,> is also symmetric, since p* and its transpose 
both satisfy (l)-(3). We should note that only p* is unique, not the 
factors Ei and Fj in (3). For instance, replacing them by EJa and aFj 
for all i,j, and any a > will produce the same p*. The extent of this 
nonuniqueness is characterized later. 
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Kruithof pointed out two desirable properties possessed by the 
projection scheme. The first, called "reversibility," says that after 
projecting traffic from q at time 1 to p at time 2, we can turn around 
and project backward to time 1, recovering q exactly. The second 
property, "divisibility," says projecting from q to p and then from p to 
r at time 3 yields just the same result as projecting directly from q to 
r. Thus the projection scheme is not noisy, in the sense of irretrievably 
losing information along the way about the initial traffic. Rather, p 
depends only on q and the row and column sums, not on the path 
followed over time. Kruithof also described a third desirable property, 
"separability," which did not hold for his basic method. The idea is to 
be able to merge or split a collection of points i orj, without affecting 
the projected traffic for any other points. By careful generalization of 
Kruithof 's method, a property similar to this can be introduced. 

III. NETWORK FLOW CONSIDERATIONS 

An important aspect of Kruithof 's method may be visualized by 
means of a flow on the simple directed graph in Fig. 1. The nodes i = 
1, 2 • • • M and j = 1, 2 • • • N represent originating and terminating 
points. The edge joining node i to nodey carries traffic p y . Interpreting 
(l)-(2) as conservation laws, the remaining edges to source s and sink 
t carry total traffic quantities 6, and dj, while the net flow from s to t 
is the sum p of all the flows p y . When an element qy vanishes in q, the 
corresponding edge from i to j is deleted in the network, so that flow 
ptj is automatically zero. For the flow p to have the form (3), all 
remaining edges must have nonzero flow p y on them. Conversely, if a 
flowpy can be constructed that satisfies (l)-(2) and does not vanish on 
any edge of the network, then it fulfills the existence conditions, so 
that (l)-(3) have a solution p*. 

The labeling method of Ford and Fulkerson 8 immediately springs to 
mind as a means of constructing a flowp v - to satisfy (l)-(2). This is a 
simple, efficient, easily programmed algorithm that maximizes the net 
flow from s to t We just assign maximum capacities of 6, for edges 
from s to i, infinity for edges from i toy, and dj for edges from; to t. If 
the maximum flow obtained is less than: 




Fig. 1 — Network model. 
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p = I bi : - I d jt (6) 



; 



then no such solution of (l)-(2) exists. The only trouble is that the 
flow obtained may not be positive on all edges, as required to establish 
that (l)-(3) have a solution. This can be cured by assigning a suffi- 
ciently small lower bound p £ > ss h > to flow on the edges from i to / 
An initial flow of h on each such edge will then be feasible. 

In practice, the total traffic quantities 6, and dj are integers or can 
be scaled up and rounded off to an approximate integer problem 
without loss of credibility. The labeling method may then be modified 
to carry h as an infinitesimal quantity, while preserving pure integer 
arithmetic for efficiency and to avoid roundoff error. Indeed, since one 
wants to draw a yes-or-no conclusion about existence of a solution, 
roundoff introduces an unwelcome uncertainty. In the modified 
scheme, all capacities and flow variables are carried as pairs (m,n) of 
integers, representing the expression m + nh. In the steps of the basic 
labeling algorithm, 8 one adds and subtracts various capacities and 
flows, attaching labels to nodes according to whether or not the result 
exceeds zero. In the modification, the obvious rules apply for adding 
and subtracting quantities m + nh; it remains, however, to specify an 
interpretation for inequalities involving such quantities. Two classes 
of inequalities must be defined: 

(i) m + nh> 0(1) means m & 1. 
(ii) m + nh> 0(h) means m > 1 or m = but n>l. 

In designating h an infinitesimal, we are really promising to choose it 
as small as necessary so that | nh \ < 1 for all n that arise. 

The solution proceeds in two phases. First the standard labeling 
method is used, but with all inequality tests to be taken as > 0(1). 
Thus, in each iteration, an augmenting path of capacity > 0(1) is 
found. The net flow increases by at least one unit (plus or minus some 
nh), so that this phase terminates after a finite number of iterations. 
At this point, the net flow must be p — nh for some h 2* 0, or else no 
flow solution of (l)-(2) exists; in effect, we have the maximum flow for 
the special case h = 0. The second phase repeats the standard labeling 
method, but with all inequality tests to be taken as > 0(h). Now each 
augmenting path has capacity > 0(h), increasing net flow by at least h 
units and, again, the iterations terminate in a finite number of steps. 
If the maximum flow reaches p, then a solution p of (l)-(2) has been 
constructed that fulfills the existence conditions; otherwise, no such 
solution exists. The two-phase algorithm was programmed, for the 
network associated with an arbitrary M x N matrix q, in about eighty 
lines of Fortran. In tests, it was able to settle the question of existence 
of solutions to specific Kruithof problems with gratifying rapidity. The 
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same technique can be used to maximize flow on any network under 
any mixture of > and >■ capacities. 

Still more can be gleaned from the network model above. Observe 
that in (6) the sum of the 6, must equal the sum of the dj in order to 
conserve flow. This requirement is just a simple example of a class of 
necessary conditions arising from (l)-(3). More generally, let /be any 
subset of originating nodes i and J{I) be the subset of terminating 
nodes j with edges to nodes in /. That is, node j is in J(I) if g,> is 
positive for some node i in I. Then the flow bi to node i in / can only 
pass to nodes in J(I), so that it is included in the flows dj out of these 
nodes. Thus the total flow Yd) from nodes in J(I) cannot be less than 
the total flow X{I) to all the nodes in 7: 

X(I) = £&,*£ Y(I) = £ dj (7) 

for each proper subset I (that is, I not empty and not containing all i). 
Conversely, these necessary conditions guarantee that every cut has a 
capacity of at least p, so that (6)-(7) are also sufficient conditions for 
existence of a flowpg that satisfies (l)-(2). 

Accounting for eq. (3), the requirement of positive flow on each edge 
allows conditions (7) to be strengthened. Indeed, flow Y(I) from «/(/) 
includes the flow X(I) into / plus any additional flow to J(I) from 
nodes i that are not in /. If any edge joins a node outside I to J(I), 
then its flow is positive and (7) becomes a strict inequality: 

X(I) < Y(I). (8) 

If equality holds in (7), the rest of the network has no connection to J 
and </(/), except through s and t Such a disconnected situation 
represents two or more independent Kruithof problems, which ought 
to be treated separately from the outset. Indeed, the rows and columns 
of qij can then be renumbered so that it is partitioned into two or more 
uncoupled blocks. To simplify the statement of later results, we assume 
that the problem does not decompose in this way, so that the network 
is connected and (8) holds for every proper subset / of the nodes i. In 
complete analogy to the case of (6)- (7), necessary conditions (6)- (8) 
are also sufficient for existence of a positive flow satisfying (l)-(2), and 
hence for existence of p*. To see this, we first reduce all network flows 
and capacities by the initial feasible flow h used in the modified 
labeling method. For sufficiently small h > 0, conditions (7)-(8) now 
apply to the reduced bi and dj as well. But these again show enough 
capacity on every cut that a flow/?,, - h 5* exists satisfying (l)-(2). 
Conditions (6) -(8) can also be used in a direct proof that p* exists, 
independently of results in the appendix. The idea is to seek a station- 
ary point of the quotient Num/Den of two multinomials in the growth 
factors Fj 2* 0. The numerator is positive on the interior and vanishes 
at boundaries (that is, where some Fj = 0): 
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Num = n Ff, 0) 

j 

while the denominator is also positive on the interior: 



Den = Y[ 



-\>>, 



2 QyFj 



i 



(10) 



Setting a derivative of Num/ Den with respect to Fj to zero yields the 
same result as substituting (4) into (5). Thus, if the values Fj make 
Num/Den stationary, they will also yield a solution p* of (l)-(3). Note 
that Num and Den are both homogeneous of order p from (6), so that 
their quotient is positive and constant along interior rays (that is, 
along aFj for all a > 0). Since the rays form a compact set and Num/ 
Den is positive and continuously differentiable on the interior, it is 
enough to show that the quotient goes to zero on the boundary to 
deduce that it achieves a (stationary) interior maximum. Now suppose 
that Den vanishes at some boundary point F and let I be the set of 
nodes i for which the factor £ QijFj in (10) is zero. Then Fj must vanish 
if qij is positive for some i in J, and hence for every j in «/(!). If Z 
measures distance from an interior point F' to boundary point F, then 
the numerator will vanish as F' approaches F at least as fast as Z Y(I) , 
while the denominator goes like Z X{I) , from (7). Thus Num/Den 
approaches zero at each boundary point F, from (8), completing the 
proof. The enterprising reader may find it instructive to ferret out the 
connection between the preceding proof of existence and the more 
general proofs in the appendix. 

One immediate consequence is that the Kruithof problem always 
has a solution p* if q is strictly positive and (6) holds. Indeed, each 
node i has an edge to every node j, so that Y(I) = p for every proper 
I and (8) follows. To verify this case more directly, one can construct 
the positive flow p,> = bidj/p, which satisfies (l)-(2). Another useful 
example is a square matrix qy with zeros only on the diagonal. Then 
Y(I) = p for every I with two or more nodes, so that only the unit sets 
i must be checked. Conditions (8) now reduce to the requirement 6, 
+ di<p for every node i. An important application involves choosing 
all bi and dj equal to one, so that M = p — N and q must be reduced 
to a doubly stochastic matrix p*. But (7) says that X(I) and Y(I) are 
just the sizes |/| and \J(I)\ of 7 and J(I), respectively. A result from 
matching theory, the "marriage theorem," 9 now shows that the con- 
ditions |/| s£ \J(I) | from (7) are equivalent to q having some positive 
principal diagonal, and that |/| < \J(I)\ from (8) imply that each 
positive element of q lies on a positive principal diagonal. These are 
the conditions cited by Sinkhorn and Knopp 5 for the doubly stochastic 
case. In a typical case, the modified labeling method will ordinarily be 
easier to apply than (6)-(8) in testing for existence of a solution. 

The growth factors E and F are essentially unique, except for the 

KRUITHOF'S PROJECTION METHOD 523 



possibility of scaling them along rays as described earlier. Indeed, let 
E and F satisfy py = qgEiFj = qijEiFj for all i and j. Then %/Ei = 
Fj/Pj = a whenever qtj is not zero, and the factor a is independent of 
I andy since the network is connected. When the problem decomposes 
into independent Kruithof subproblems, the network breaks into sep- 
arate connected components, each with a separate scaling factor a. 
The case that p and q are symmetric implies that E = a¥, so that the 
growth factors may be scaled to satisfy E = F uniquely. 

IV. CONVERGENCE CONSIDERATIONS 

In iterating (4) and (5) to solve for p*, a good convergence test can 
be based on the norm consisting of the sum of the absolute values of 
the differences between the left and right sides of (1) and (2): 

g - S I bi - E t I gift | + I \dj - Fj £ qijEi | . (11) 

«' j j « 

Clearly, g is the net error in traffic units and goes to zero as p goes to 
the solution p*; we can show that, in fact, g does so monotonically. 
Indeed, g is continuous, convex, and piecewise linear in E, with Ei- 
derivatives of the form: 

- I QuFj 

J 



sgn bi - Ei £ QifFf + sgn dj - Fj X q f jBi 



and this expression always takes the sign of Ei £ qijFj — bi or else 
vanishes. Thus g can only decrease or remain constant as Ei is 
increased or decreased to satisfy (4) and g achieves its minimum, for 
any fixed value of F, when each Ei is given by (4). A similar discussion 
holds with respect to the F y -derivative of g. Accordingly, when g 
becomes small during iteration, it will remain so, and it is appropriate 
to stop. 

During the process of iterating to compute E and F, we can accu- 
mulate a value of g with very little additional effort. Specifically, when 
all the Fj have a new value, the second term of (11) vanishes. Proceed- 
ing to update the Ei, we compute X QijFj for use in (4). With two more 
additions and one multiplication, we get the corresponding contribu- 
tion to g in the first term of (11). By the time all new E t are computed, 
a value of g for F and the old E is available. Clearly, the iteration may 
be interpreted as a relaxation scheme to minimize g by cyclically 
minimizing over the Ei and Fj. 

Since nonexistence of a solution p* when (6) holds can only accom- 
pany zero elements in q, one might attempt to force a solution by 
substituting a small positive value for each zero. In general, this is a 
terrible idea; Sinkhorn, 4 for instance, shows some examples of patho- 
logical behavior associated with such schemes. The iteration process 
will seek to get significant flow on some edges with small qy by using 
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very large growth factors. The solution p* will not greatly resemble q 
and the iteration will ordinarily converge very slowly. Indeed, slow 
convergence is a possible warning that the problem is ill-conditioned 
in some way. For such cases, it is prudent to check the "existence 
margin" by setting successively larger elements of q to zero and 
running the existence test until it fails — just the reverse of adding 
small positive terms. The labeling method is sufficiently economical of 
computing time that it may be repeated often in preference to perform- 
ing a great many iterations. 

In general, ill-conditioning occurs when some of the total traffic 
values bi and d, are not particularly consistent with one another. That 
is, when some collection of small elements qtj are set to zero, (l)-(2) no 
longer have any solution for which p# vanishes if and only if q$ does. 
The limit as these elements approach zero may not exist, or it may 
depend on the specific way in which they vanish. In short, p* is not 
necessarily continuous in the elements of q at the boundaries of the 
feasible region. On the other hand, p* is provably continuous in the 
row and column sums, or any other constraint levels, so that adjusting 
them to achieve consistency at the boundary is a stable procedure. 
Thus, the proper way to treat ill-conditioning is by readjusting some 
values of 6, and dj to make them more consistent, though it may not 
be obvious how to do this. We see later that there is an easy way to 
extend the Kruithof method so that it automatically allocates traffic 
among the rows or columns of prespecified aggregations in a consistent 
and reasonable way. 

It is also possible for a quite reasonable problem to converge at a 
very slow rate. For example, a problem may decompose into multiple 
independent subproblems, whose network components are only con- 
nected by way of s and t. Now, each of the subproblems may be well- 
behaved, so that the overall iteration process converges rapidly. Nev- 
ertheless, when a few small positive values of qtj are introduced to 
couple the subproblems, convergence will be rapid at first and then 
become rather slow, as a rule. Moreover, the solution to which the 
problem converges is a sensible one that differs only slightly from the 
decoupled case. This model could apply to two or more countries, for 
example, with much more traffic internally than across their borders. 

What causes the above difficulty is the difference in degree of 
uniqueness between coupled and decoupled cases. For n subproblems, 
n independent arbitrary scaling factors a will appear in the general 
solution. The actual values they assume will be determined by the 
initial values assigned to E or F. This effect appears as some arbitrar- 
iness in the relative sizes of the E, for those portions of E associated 
with the various subproblems. For the coupled case, only a single 
overall scaling factor is appropriate; the portions of E from the different 
subproblems must now be scaled in a correct ratio to each other. In 
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the rapidly convergent phase, each subproblem is solved to within its 
scale factor; the slow phase corresponds to a process of adjusting scale 
factors to account for small coupling terms. 

The slowly convergent phase may be shortened appreciably by 
means of standard numerical schemes for acceleration of convergence. 
Wynn's algorithm, 10 in particular, has been employed successfully to 
project the values of the Ei for successive steps, in order to estimate 
their limit. The strategy is to calculate the norm g for the projected E 
and, when this is much less than the current error g (one-fiftieth, say), 
restart the iteration from the projected value. Additional computa- 
tional effort and program steps for a convergence acceleration option 
in Kruithof 's method are minor. Of course, acceleration techniques 
cannot rescue a truly ill-conditioned case, where the existence margin 
is small. All this is not meant to imply that slow convergence is the 
rule with the Kruithof method. In fact, some rather large examples, 
involving several hundred rows and columns, have been solved quite 
readily. 

V. VARIOUS EXTENSIONS 

An immediate generalization of the basic Kruithof scheme would be 
to stratify the traffic data in more than two dimensions. Besides the 
originating and terminating points i and j, other indices k, I, m 
might specify time of day, week, or year, type of traffic (business or 
residential, for instance), and so on. A three-dimensional case of the 
general Kruithof problem then might take the specific form: 

p ijk = q ijk EiFjG k (12) 

b i = Y iP i jk = E i Y 1 qijkFjG k (13) 

dj = £ Pijk = Fj J qu k EiGk (14) 

i,k i,k 

fk = %Pijk = G k 2 qijkEiFj. (15) 

i.j 'J 

Reduction to the standard case in the appendix proceeds as before. 
We define events e = ( i, j, k ) and constraints c corresponding to each 
value of i, j, and k, while p and q are normalized to probabilities by 
dividing by the sums of their elements. 

The proofs of existence, uniqueness, and convergence in the appen- 
dix still hold for this case. In particular, a solution p* of (12)-(15) exists 
if and only if (13)-(15) have a solution />,>■* that is positive or zero 
accordingly as g#* is positive or zero. The norm g consisting of the sum 
of absolute values of differences between left and right sides in 
(13)-(15) is still net error in traffic units, though it no longer decreases 
monotonically with each new value of Ei, Fj, or G k . Nevertheless, g 
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goes to zero as p goes to p* and it is still a reasonable indicator of 
convergence. When q is strictly positive, the interior solution />,,* = 
bjdjfk/p 2 of (13)-(15) demonstrates that p* also exists, provided that 
the bi, dj, and /* each add up to p. A more general existence test would 
involve solving the linear programming problem described in the 
appendix, with appropriate precautions against being misled by the 
effects of any roundoff error. Network flow models no longer apply, 
and integer solutions need not occur in the case of integer constraints. 
Another three-dimensional example of the general problem might 
be as follows: 

Pijk = qijkEijFjkGik (16) 

bu = S p^ = Eij £ qijkFjkGik (17) 

* * 

djk = £ Pijk = F Jk £ qijkEijdk (18) 

i i 

fik = I Pijk = Gik £ qukEijFjk, (19) 

/ j 

and the same sort of discussion applies to this case as to (12)-(15). For 
four or more dimensions, the reader should have no difficulty creating 
a great many extensions of this general class, such as pyu = 
QijkiBijCjkDkiEuFikGji. The rule is to multiply (fa*/... by one factor 
Eij. . . for each constraint in which />,>*/. appears, and then solve each 
constraint for its factor by dividing into the constraint level 6,> . In 
higher dimensions, the number of elements in p and q grows much 
faster than the numbers of constraints and multipliers. The numbers 
of the latter thus remain reasonable, if only to stay within storage and 
computational limits on the former. The corresponding linear program 
to test for existence will therefore have a basis of reasonable size, as 
well. 

Another class of extensions involves specifying less about total traffic 
quantities in the two-dimensional case (or any other dimension, using 
the previous generalization). That is, various sets I or J of originating 
or terminating points i or j may be lumped together, with only their 
total traffic D and d to be supplied externally. Now (l)-(2) yield the 
following constraints: 

b, = Zb, =ZZpu (20) 

/ / j 

dj = %dj='Z'Zp lJ . (21) 

j j i 

Results on the general Kruithof problem in the appendix show that Ei 
= Ei and Fj = Fj in this case, so that (3) becomes: 

Pij = q.jE.Fj = qijE,Fj, (22) 
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assuming each i belongs to one J and each j to one J. 

To solve for £ and f , we can collapse the problem to a simpler form: 
add together all rows i in J and all columns j in J to work with reduced 
matrices p and q, as follows: 

pu^llPij ffi'-SZtftf- (23) 

i j i j 

Now (20)-(22) are reduced to exactly the same form as (l)-(3): 

5/=EP/</ dj^^pu pu=quEiPj, (24) 

j i 

but with many fewer constraints to be satisfied. Note that this proce- 
dure of aggregating traffic nodes may absorb some elements of q that 
were zero or small, in order to ameliorate Ul-conditioning. Indeed, the 
corresponding disaggregation formula (22), to be used after we have 
solved (24), simply allocates traffic py to the element qy in proportion 
to its relative contribution to qjj in (23). This in turn automatically 
shares out bi and dj among their 6, and dj in a consistent manner, as 
mentioned earlier. As another of its virtues, aggregation is a smoothing 
process that can cut down the effects of errors in the predictions of 
total traffic by reducing the number of independent parameters. The 
reduction in manual and computational effort is also an evident 
advantage. To organize the computation in an efficient manner, we 
would start with two tables, I(i) and J(j), that assign each point i or 
j to its appropriate aggregate. Then we run through the pairs (i, j), 
adding each qy into its correct qmjy\. After we solve (24) for £ and 
f\ the answer is just/ty = qyEm&JU) from < 22 )- Existence testing can 
be performed directly on q, b, and d with the labeling method. 

The scheme above illustrates a sense in which a "separability" 
property can be introduced, similar to what Kruithof sought. There is 
no real need to require that different / or J be disjoint. If overlap is 
permitted, then qy would be multiplied by Ei for each J that contains 
i, and by Fj for each J containing/ Collections of rows or columns can 
now be aggregated only if they all he in the same sets I or J. Since 
some pa may now be counted twice, the constraints can no longer be 
interpreted as conservation laws for a flow, and the existence test 
becomes a linear program, as described in the appendix. 

Another way of specifying less in the Kruithof problem is to leave 
some rows or columns unconstrained. At such i and/', we can specify 
the growth factors Ei and Fj arbitrarily; Bear 2 has considered choosing 
these multipliers to be one. One advantage of not constraining is that 
any element qy for which row i and column/ are not constrained plays 
no part in the solution process and may be set to zero for convenience. 
For computational efficiency, we multiply each of these rows by its 
fixed growth factor and add together all such rows to form a single 
new row; similarly, all unconstrained columns combine. Of course, 
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fewer specifications may be introduced in higher-dimensional schemes 
as well. 

Anyone can tailor their own ad hoc constraints to account for 
additional knowledge of the future. For instance, it may be known that 
certain items of point-to-point traffic py are growing considerably 
faster than the other, more typical items in their rows and columns. 
Then a constraint may be created to fix the projected value of the sum 
of all such items. This produces one more growth factor to be multi- 
plied into qij for these selected items only. Similar treatment may be 
given to a class of items that have slower than normal growth, that 
decrease, or that just vanish. The general Kruithof problem defined in 
the appendix includes all such cases, as well as any other linear 
constraints that may need to be introduced. 

VI. SUMMARY AND COMMENTS 

In this study, we have sought to give a unified view of Kruithof 's 
teletraffic projection method, including theoretical aspects as well as 
practical details for its implementation. The mathematics in the ap- 
pendix treats the theory of a general Kruithof problem. Necessary and 
sufficient conditions for existence and convergence of its solution are 
derived, along with proofs of uniqueness and continuity. The main text 
considers special cases of this problem that are of particular interest. 
Schemes for existence and convergence testing and for handling slow 
convergence are discussed. We conclude by trying to place this projec- 
tion scheme in the context of the Bell System planning function. 

An early and important step in the Bell System planning process is 
that of predicting future demand for the various services offered. By 
their nature, such projections can be quite uncertain, since they will 
include cumulative effects of several years' fluctuations in the United 
States and world economy, for instance. Indeed, analysis of time series 
of typical traffic data" indicates that about five percent per year of 
random error remains in even the best projections, and must be 
regarded as inherently unpredictable. Nevertheless, a strategy is avail- 
able to cope with such uncertainties, for the purposes of planning. 

The fundamental assumption required in this strategy is that traffic 
increases monotonically with time. First a plan can be generated, based 
upon some "best guess" of the demand profile over time. From year to 
year, the time scale of the plan can then be corrected to match up the 
originally projected demand with actual values or better estimates, 
based upon more recent data. In effect, a parameter such as total 
traffic is thus used as a new independent variable in the plan, while 
time is a dependent variable that absorbs much of the economic 
fluctuations and other error. However, this leads us to view the overall 
process of planning as a system, rather than a collection of independent 
modules, one of which is projection. We see that a "sliding time scale" 
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approach to planning now places a premium not so much on absolute 
accuracy of traffic predictions over time as on "relative accuracy" or 
consistency and uniformity among the various traffic quantities that 
are projected. 

A typical plan may be based upon many thousands or tens of 
thousands of projected traffic items. We wish to predict these quanti- 
ties such that a single readjustment of the plan time scale (based on 
some total traffic measure, for instance) can do a reasonable job of 
correcting for the error in each item. This requires that projection be 
done in such a way that the prediction errors in individual traffic items 
will tend to be highly correlated. Thus, a scheme which projected 
individual time series for each separate traffic quantity, for example, 
might give the best absolute accuracy for each item, but still be 
unsuitable for planning purposes because the noise components in 
these time series would tend to be independent. At the opposite 
extreme, initial measurements of all traffic quantities might simply be 
increased by a single overall growth factor for each year under study. 
This would produce very strong correlations but fails to take into 
account detailed knowledge of growth patterns, say, for separate 
portions of a study area. 

The general Kruithof method offers many middle roads. Any collec- 
tion of average or overall traffic quantities b may be predicted exter- 
nally (from time series, for example, or market surveys). As shown in 
the appendix, the remaining items can then be projected to be con- 
sistent with whatever is given. Inserting a great many external predic- 
tions introduces more detailed knowledge, but reduces the correlation. 
Supplying fewer external specifications yields stronger correlations, at 
some loss in accuracy; various tradeoffs are possible. 

All the schemes of Kruithof type act to minimize the net information 
change in the projection, subject to those external constraints being 
enforced. This gives them the remarkable properties called "reversi- 
bility" and "divisibility" by Kruithof. Essentially, all that is lost of the 
original data q in projecting it to future values p is whatever is inherent 
in the externally provided average quantities b. Supplying new values 
b' for these quantities will thus allow us to continue the projection 
from p to another year or recover the base data q exactly. Effectively, 
Kruithof's method is able to resolve q and p into a part which 
determines some arbitrarily chosen system of average quantities b 
that are to be changed and an "orthogonal" part that does not change. 
This latter part is essentially the equivalence class C(Q) discussed in 
the appendix. 

Kruithof's original proposal, and much of the subsequent work on 
the subject, is concerned with projecting two-dimensional arrays/),; of 
traffic data. In this case, the most natural overall quantities to be 
specified externally are the sums of rows i and columns j in p. For 
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example, supplying two hundred parameters for a 100 X 100 matrix 
would suffice to project a total of ten thousand items. An option is to 
aggregate rows and columns, perhaps in collections of average size 
four, so that only fifty external sums are needed. This is a tradeoff 
that increases correlations but may decrease accuracy. Meanwhile, it 
can alleviate possible ill-conditioning, reduces manual and computa- 
tional effort, and still projects ten thousand items. Examples of data 
organized in three or more dimensions are also known in Bell System 
planning. Kruithof's method generalizes to such cases without any 
particular difficulty. 
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APPENDIX 

In this appendix, we define the general case of the Kruithof problem 
and derive various properties. In particular, necessary and sufficient 
conditions for existence and convergence of solutions are developed, 
and uniqueness and continuity are proved. To begin with, consider two 
probability distributions, P e and Q e , over the same finite set of disjoint 
events {e}, so that: 

£ P e = 1 P e 2* (25) 

e 

ZQe=l Qe^O. (26) 

e 

The information content of a probability, in appropriate units, is minus 
its logarithm. Thus the change in information from Q e to P e is just log 
P e — log Qe = log( Pe/Qe). The average of this information change is 
defined as: 

JC(P,Q)-ZPelog(P,/Q.) f (27) 

e 

which will be interpreted as a measure of how close distribution P is 
to distribution Q. A simple example is the case that all probabilities Q e 
are equal; now K reduces to a linear function of the entropy of 
distribution P. Thus entropy measures departure from the equiprob- 
able case (corresponding to classical equilibrium). The expression (27) 
goes by several names in the literature; for instance, Kullback distance, 
/-divergence, relative entropy, discrimination information, and Gibbs 
free energy. 
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As a distance measure, K has two desirable properties: it is not 
negative and vanishes only when P = Q. Two other properties, sym- 
metry and the triangle inequality, are lacking; Csiszar 7 points out, 
however, that laws analogous to the parallelogram identity and Pytha- 
goras' theorem do hold, with K playing the role of squared length. To 
show that K is nonnegative, first note that F(X) ■ log ( 1/X) is strictly 
convex, since F" = 1/X 2 > 0. Defining X, = Q e /Pe and using convexity 
in (27) yields the inequality: 



K = X PeF(Xe) 2* f( £ PeXA = F(l) = 0, 



(28) 



with equality only if all P e F(X e ) vanish, in which case P = Q from 
(25)-(26). 

When Q e is positive and P e approaches zero, P e log( P e /Qe) goes to 
a zero limit; to ensure continuity, we will define it to vanish at P e = 0, 
even for the case Q e = 0. If some Q e = for positive P e , then K is 
infinite. Confining our attention to P that are not infinitely far from Q, 
P e must vanish whenever Q e does, so that e is an event of probability 
zero. With no loss of generality, such e may be excluded from the set 
{e} of events for now, so that all Q e are positive. 

Consider the problem of minimizing K(P, Q) for fixed Q, over all 
distributions P subject to (25) and a finite set {c} of arbitrary linear 
constraints having the general form: 

ZP e A ec = B c . (29) 

e 

This amounts to finding the distribution P that is closest to Q on the 
intersection (denoted S(B), or just S) of the positive orthant P 2* and 
the hyperplanes (29), which prescribe that certain averages over P 
take on the values B c . (To simplify the notation, assume that the 
equality in (25) is designated c and is included among the constraints 
c.) Observe that S is a compact convex polytope, so that the continuous 
function K achieves its minimum value on S, whenever S is not empty. 
Further, this minimum occurs at a unique point P* in S, and there are 
no other local minima of K, since it is strictly convex. To see this, note 
that the matrix of second partial derivatives of K with respect to P is 
diagonal and positive definite. We assume from now on that S contains 
more than one point. 

Suppose that S contains an interior point P Inl (that is, no Pi"' 
vanishes); then P* is also an interior point. Indeed, consider any 
boundary point P Bdy and the line p flrf >P /n ' joining it to P /n '. The 
gradient of K has components 1 + log (P e /Qe) that become arbitrarily 
large negative on some neighborhood of P Bdy for those P Bdy that 
vanish, while all other components remain bounded. Thus, a segment 
f pBdypim containing P Bdy can be found along which K decreases 
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toward the interior, and P Bdy cannot be a minimum point of K. Since 
the gradient of K is continuous on the interior of S, an interior P * is 
a stationary point of K and, by strict convexity, there are no other 
stationary points of K on S. 

When P* is an interior point, the minimization problem can be 
solved by using stationarity. Specifically, adjoin constraints (29) to K 
with multipliers v c to form the Lagrangian function L ( P, v ) , as follows: 



L = lP e 



\Og(Pe/Qe)-ZAeeV c 



+ I BcVc. (30) 



Now L is strictly convex on the positive orthant P ^ and becomes 
infinite if any P e does, while its gradient is negative infinite on the 
boundaries. Thus it achieves a unique minimum over P 2* at some 
interior stationary point P(v) > for any fixed values of v c . This point 
is found by requiring the P e -derivative of L to vanish for each e, 
yielding the following necessary conditions: 

w e m \og(P e /Q e ) = X A ec v c - 1. (31) 

Setting V c = exp(v c ) in (31), except for the constraint c corresponding 
to (25) with Vc = exp(t;<- — 1), now yields 

P.. = Q e exp(^ t .) = Q ( . n V?-% (32) 

c 

with all V r strictly positive, so that P(v) cannot be on the boundary. 

Suppose some v c are found such that P(v) satisfies the constraints 
(29) and thus lies in S. Then since L = K on S, P(v) is a stationary 
point of K on S, and hence is the unique minimum point P*. Con- 
versely, the linear program of minimizing dK for all small variations 
dP that satisfy (29) has (31) as its dual constraints. Since the primal 
problem has rfP = as an optimum at P*, the dual is feasible there, 
sc that v can be found to satisfy (31) at P*. The point of all this 
reasoning is that a solution to (29) can be found with the specific 
product form (32) if and only if S possesses an interior. When such a 
solution exists, it is also the unique minimum point P * of K over S. 
We can now define the Kruithof problem, in general, as that of finding 
the factors V c in (32) so as to satisfy the constraints (29). 

Kruithof's "reversibility" property amounts to symmetry of the 
"closeness" relation: whenever P is closest to Q, then Q is closest to P 
in the same sense. That is, we define some new constraint levels: 

Be m I Q e A ec (33) 

and seek a distribution R to minimize: 

K = K(R, P) = £Re \og(Rc/P e ) (34) 
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over the set S(6) defined by linear constraints: 

£ R e A ec = 6 C R e & 0. (35) 

e 

Defining V c ■ 1/V C , we can write Q e from (32) in the product form: 

Q e = Pe II flK (36) 

From (26), (33), and (36), Q is an interior point of S(6) having the 
product form (32), so that R = Q is the unique minimum solution for 
R, and Q is closest to P. Kruithof's "divisibility" property is just 
transitivity of the closeness relation: when P is closest to Q and R is 
closest to P, then R is closest to Q. Specifically, we choose an arbitrary 
constraint vector fi in (35), such that S(B) has an interior point, and 
seek R to minimize it in (34) over S(6). The solution for R is the right 
side of (36) for some V, and substituting for P from (32) yields: 

Re = QeT[(V c V c ) A », (37) 

c 

which again has the product form (32). Thus R minimizes K(R, Q) 
over S(6) and is closest to Q in this sense because it is closest to P. 

Symmetry and transitivity show that "closest to" is an equivalence 
relation determined by the particular matrix [A K ]. This relation 
partitions the set of positive distributions over {e) into equivalence 
classes. Each class C(Q) can be generated from any one of its members 
Q by using (32): choose all positive values of the V c for c ^ c and scale 
the resulting values P as necessary to meet the normalization condition 
(25). Since the column of [A ec ] corresponding to c is all ones, V c 
appears linearly in (32), and the scaling above represents a particular 
choice of that variable. Uniqueness says that each constraint vector B 
is achieved at most once in each class, while the existence condition 
asserts that B is achieved in every class if it is achieved in one class. A 
natural mapping B = /"(P), namely the linear mapping (29), takes any 
C(Q) into the set E of constraint vectors B for which S(B) has an 
interior. Clearly, /"is continuous and is one-to-one and onto E from the 
existence and uniqueness results. We will see later that /is one-to-one 
on the closure of C(Q), which is compact from (25). (However the 
closure is not necessarily an equivalence class.) It follows that / is a 
homeomorphism of C(Q) and E. In particular, P* is a uniformly 
continuous function of the constraint vector Bon£ and its closure. 

A useful result follows from the linear relation between v and w in 
(31). Specifically, let R be any solution of (25) and (29), multiply R e by 
w e , and sum over e, using (29) and (31) to obtain 



£ R eWe = %Re 



^AecVc - 1 



= ZB c v c -l (38) 
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by exchanging the order of summation. It is easy to show from (31) 
that the left side of (38) is just K(R, Q) - K(R, P), while the right side 
does not depend on the particular choice of R in S(B), so that 

K(R, Q) + K(P, P) = K(R, P) + K(P, Q) 

for all R, P in S and all P, Q in C. But now the choice P = P* - P 
yields an example of the Pythagorean theorem noted by Csiszar: 

K(R, Q) = K(R, P*) + K(P*, Q). (39) 

Roughly, it says that the distance from R to Q breaks into a component 
within S from R to C(Q) at P* and an "orthogonal" component from 
P* to Q within C(Q). 

Now we will investigate the nonlinear (Wolfe) dual problem to 
minimization of K on S. Let H (v) be the minimum of L(P, v) over the 
positive orthant P > 0, so that 

H(v) = L(P(v), v) *£ L(P*, v) = K(P*, Q) < K(R, Q), (40) 

which shows that H (v) is bounded above if S is not empty. Substituting 
(31), (32), and (38) into (30) allows us to express H as a function of 
each of v, V, w and P as follows: 



H = I BcV c -%Qe expl £ AecV c - 1 

c e \ c t 

= 1 + I B c l0g( V c ) - I Q e n Vi" (41) 

c e c 

= 1 + £ [R e We ~ Qe exp(lV e )] = 1 + % [Re log(P e /Qe) ~ Pe]- 

e e 

Direct differentiation of H(P) shows that it achieves a unique maxi- 
mum over the positive orthant at P = R where H = K(R, Q). (Break 
H into a linear part for those components of R that vanish and a 
strictly concave part.) Indeed, if H(P) goes to K(R, Q) on P 5= 0, we 
can conclude that P approaches R. Now the difference between K(R, 
Q) and H (P) achieves a unique minimum of zero at P = R: 

K - H = X [Re log(R e /Pe) + Pe] - 1 = K(R, P) + £ P e - 1. (42) 

e e 

Whenever this expression vanishes for some P of the form (32), we 
have K(P*, Q) = K(R, Q) from (40), and thus P* = R by uniqueness 
of the minimum. 

Suppose that values of P(v) of the form (32) approach some limit R 
on the closure of C(Q). Such R satisfies (25) and thus (29) for the 
constraint vector B = f(R). Now K - if becomes arbitrarily small, and 
we conclude that R is the P * that minimizes K. This shows that the 
closure of C(Q) consists of points P* that minimize K(P, Q). The 
uniqueness of the minimum then says that the mapping /is one-to-one 
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on the closure, as promised earlier. Choosing R = P* in (41) shows 
that H(P*) = K(P*, Q) even when H(P) is extended to the closure of 
C(Q) by continuity. Finally, the result (39) again follows by setting P 
= P*in(42). 

The nonlinear dual problem consists of maximizing if over all v, and 
thus over positive V, or over the linear affine space of values w 
generated by (31). By its construction, H is concave in v, while 
differentiation shows that it is strictly concave in each individual v c . 
Setting the Vc-derivatives of H to zero yields necessary conditions for 
a stationary maximum: 

B c = £ QeA ec n Vfr\ (43) 

e c' 

These are the same relations that would be obtained by substituting 
(32) into (29), namely, the general Kruithof problem. One possible 
scheme for solving eqs. (43) would be by relaxation: choose some 
variable V c and adjust it to maximize H with all other variables V c - 
held fixed; then choose some other variable and repeat, cycling through 
all V c infinitely often. (When H is bounded above by K(R, Q) in (40), 
a value of V c to maximize H and satisfy (43) always exists uniquely, 
because —H is strictly convex in each v c and is arbitrarily large for 
large v c .) If the values of P(v) for the iterates v approach the limiting 
value P*, then the relaxation procedure represents a means of com- 
puting P*. More generally, any collection of constraints c could be 
solved simultaneously in (43), followed by another collection, and so 
on, so that each c appears infinitely often. Simultaneous solution of 
several constraints can be harder than the scheme of treating one 
variable at a time, however. Another possibility might be to solve (43) 
approximately for V c , so that H increases at each step, but not 
necessarily to its exact maximum in V c . The general relaxation pro- 
cedure is just an attempt to maximize H over all variables by doing a 
few at a time. Such relaxation schemes are known 12 to converge to the 
maximum of a concave function under very general conditions. 

In practice, A ec will generally be a zero-one matrix, so that the 
powers of V c in (32) do not become a nuisance. In such a case, the 
constraints (29) have a simple interpretation, since they assign prob- 
ability Bato the event c that is the disjoint union of those e for which 
A ec = 1. Pursuing this view, the sum in (29) is taken over all events e 
included in c (denoted e C c) and the product in (32) is taken over all 
events c that include event e (denoted c D e). Substituting (32) into 
(29) now yields 

Be = X Q e II Vc =VcZQe II Vc, (44) 

etzc c'Oe eCc Cf*c'De 

as the form taken by (43) in this case. The relaxation iteration step 
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now reduces to dividing B c by the sum on the right in (44), in order to 
calculate the new value of V c that maximizes H. Under certain weak 
restrictions on the sequence of variables chosen for iteration, it can be 
shown that H increases to the limit K(P*, Q) and that P(\) converges 
to P*. Specifically, assume that the iteration scheme includes an 
infinity of intervals of length M , for some sufficiently large M , in each 
of which H is maximized over every variable V c at least once. Then 
the relaxation iteration converges if and only if S is nonempty; the 
limit is P*, which minimizes K over S and lies on the closure of C(Q). 
Indeed, with S nonempty, the consecutive values of H are nonde- 
creasing but bounded above in (40). Thus H approaches some limit 
H*, and the successive increases rfHmust eventually go to zero. Direct 
computation from (41) and (44) now yields the relation: 

dH = B c [dv c - 1 + exp( - dv c ) ] , (45) 

where dv c is the corresponding change in v c . Differentiation shows this 
expression to be strictly convex in dv c (except in the trivial case B c = 
0), vanishing only at its minimum, namely at dv c = 0. Thus each du c 
also goes to zero as H approaches H*, so that dw goes to zero from 
(31). Consider the values of P(v) that are obtained each time the 
constraint c corresponding to (25) is satisfied in one of the postulated 
intervals of length M . Since these values are confined to the simplex 
defined by (25), they have a subsequence that converges to some limit 
point R. Now each member of the subsequence differs from a solution 
of constraint c by at most M changes, each of order dw. It follows that 
the limit R will satisfy every constraint c . But then, from the discussion 
after (42), H goes to K(R, Q) on the subsequence and R is P*. Finally, 
H increases to its maximum over P 2= for all iterates P, so that they 
converge to P*. Conversely, whenever the P converge, the limit R 
satisfies all constraints, so that S is nonempty. Csiszar proves conver- 
gence for cylic iteration on collections of constraints, if each collection 
contains c, by an elegant application of (39). In general, such cases do 
not include single-variable relaxation schemes, such as those treated 
above. 

The artificial restriction that no Q e may vanish can now be dropped, 
since (32) shows that P e is zero whenever Q e is in any case. The 
definition of S and its interior must be modified to account for all such 
conditions, of course. That is, (25) and (29) are supplemented by 
requirements that P e vanish whenever Q e does, while P e = makes P 
a boundary point of S only if Q e is positive. Thus the Kruithof solution 
exists if and only if some P satisfies the constraints and vanishes for 
exactly the same events that Q does. 

The question of whether the existence condition is met for a partic- 
ular constraint vector B can be resolved, in principle, by constructing 
such an interior P with standard linear programming techniques. 
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Indeed, consider the problem of maximizing A?0 subject to linear 
constraints: 

£ (T e + h)Aec + u c = B c TeSsO u c 2= 0, (46) 

where each equality c has been written such that B c > 0. The slack 
variables form the initial basis u c = B c , and a phase one procedure 
minimizes their sum. If not all u c are forced to zero at optimum, then 
S(B) is empty. Otherwise, a point in S has been constructed, so that 
the relaxation iteration will converge. In this case, all the u c are 
dropped and phase two proceeds to maximize h. As soon as some step 
causes h to exceed zero, P e = T e + h is the desired interior point. If the 
optimum still has h = 0, then S(B) has no interior, and the iterative 
solution will not take the product form (32). 
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